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Quantum  mechanics  allows  for  many-particle  wavefunctions  that 
cannot  be  factorized  into  a  product  of  single-particle  wave- 
functions,  even  when  the  constituent  particles  are  entirely  distinct. 
Such  ‘entangled’  states  explicitly  demonstrate  the  non-local  char¬ 
acter  of  quantum  theory1,  having  potential  applications  in  high- 
precision  spectroscopy2,  quantum  communication,  cryptography 
and  computation3.  In  general,  the  more  particles  that  can  be 
entangled,  the  more  clearly  nonclassical  effects  are  exhibited4,5 — 
and  the  more  useful  the  states  are  for  quantum  applications.  Here 
we  implement  a  recently  proposed  entanglement  technique6  to 
generate  entangled  states  of  two  and  four  trapped  ions.  Coupling 
between  the  ions  is  provided  through  their  collective  motional 
degrees  of  freedom,  but  actual  motional  excitation  is  minimized. 
Entanglement  is  achieved  using  a  single  laser  pulse,  and  the 
method  can  in  principle  be  applied  to  any  number  of  ions. 

Most  experimental  demonstrations  of  entanglement  to  date  have 
relied  on  the  selection  of  data  from  random  processes,  such  as  the 
preparation  and  detection  of  photon  pairs  in  parametric  down- 
conversion7-9  or  of  atoms  in  a  thermal  beam10.  All  methods  of  this 
type  suffer  from  inescapable  signal  degradation  when  entanglement 
of  larger  numbers  of  particles  is  attempted,  as  the  probability  of 
randomly  generating  the  appropriate  conditions  decreases  expo¬ 
nentially.  For  instance,  in  the  experiment  of  ref.  7,  two-photon 
entangled  states  could  be  generated  and  detected  at  a  rate  of  roughly 
1,000  per  second,  three-photon  states  at  a  rate  of  30  per  hour,  and 
four-photon  states  at  an  extrapolated  rate  of  several  per  year. 
Trapped  ions  have  been  suggested  as  a  system  in  which  such  effects 
might  be  avoided11,  and  we  have  demonstrated  two-particle  entan¬ 
glement  in  a  deterministic  way12.  By  ‘deterministic’  we  mean  that 
the  desired  state  could  be  produced  with  a  high  degree  of  certainty  at 
a  user-specified  time13,  which  is  necessary  for  avoiding  the  degrada¬ 
tion  described  above.  However,  that  experiment  relied  on  the 
particular  behaviour  of  two  ions  in  a  quadrupole  radio-frequency 
trap,  and  could  not  easily  be  applied  to  larger  numbers  of  particles. 

The  entanglement  technique  proposed  by  Molmer  and 
Sorensen6,14  can  be  understood  by  considering  a  pair  of  spin-half 
charged  particles  confined  together  in  a  harmonic  potential.  The 
energy  levels  of  this  system  are  illustrated  in  Fig.  1,  where  fico0  is  the 
internal  energy  splitting  of  each  particle,  and  v  is  the  oscillation 
frequency  of  a  particular  collective  mode  of  the  particles  in  the  trap. 
Using  laser-cooling  and  optical-pumping  techniques15,  the  particles 
are  initially  prepared  in  their  spin-down  internal  state  and  in  the 
ground  state  of  their  collective  motion:  \if)  =  |  II  0).  By  applying 
optical  fields  oscillating  at  co0  +  v  —  8  and  oj0  —  v  +  5,  the  two-step 
transition  from  |||  0)  to  |||  0)  is  driven.  For  sufficiently  large  5,  the 
intermediate  states  |||  l)  and  |||  l)  are  negligibly  occupied,  so  that  no 
motional  excitation  occurs.  The  resulting  interaction  hamiltonian, 
in  the  rotating-wave  approximation  and  the  Lamb-Dicke  limit,  is 
then 


H  =  ^(  irrxill  +  III  X  TTI  +  ITIXITI  +  lltxtll)  (1) 

with  Cl  =  i]2 Cl2lb  when  the  single  particle  || )  <-*■  || )  transition  has 


Rabi  frequency  Cl  and  the  Lamb-Dicke  parameter  is  rj.  For  an 
excitation  involving  momentum  transfer  hk  and  a  total  particle 
mass  of  M,  r)  is  given  by  (hk1  IlMv)112 .  Entanglement  is  achieved  by 
applying  H  for  a  time  t  =  idlCl,  making  the  spin  wavefunction 
\\f2)  =  (HI )  —  z|||  ))/(2.  This  spin  state  is  in  fact  created  for  any 
initial  motional  state  |  n),  so  long  as  the  Lamb-Dicke  criterion 
r/2n  <C  1  is  satisfied. 

In  order  for  the  intermediate  states  |||  l)  and  |||  l)  to  be  negligibly 
occupied,  the  detuning  5  must  be  large  compared  to  the  transition 
linewidth  qCl .  However,  it  is  clear  from  the  expression  for  Cl  that  the 
entanglement  speed  is  maximized  for  small  5,  and  in  fact  the 
technique  can  still  be  applied  for  8  ~  qCl  (refs  14,  16).  Although 
motional  excitation  does  then  occur  to  some  degree,  for  select 
values  of  8  the  excitation  vanishes  at  precisely  the  time  that  the 
entangled  spin  state  is  created.  The  condition  for  this  to  occur  is 

blqCl  =  2\fut  (2) 

for  any  integer  m,  and  the  maximum  excitation  during  the  pulse 
then  has  mean  quantum  number  h  =  1/2 m.  Our  experiment  is 
operated  with  m  =  1 . 

As  discussed  in  ref.  6,  the  entanglement  method  is  scalable  in  the 
sense  that  precisely  the  same  operation  can  be  used  to  generate  the 
N-particle  entangled  state 

fc>  =  (ll  |...T)  +  iN+1||  | ... | (3) 

if  N  is  any  even  number,  while  for  N  odd,  |  \pN)  can  be  generated  using 
one  entanglement  pulse  accompanied  by  a  separate  independent 
rotation  of  each  particle’s  spin. 

If  the  ions  are  uniformly  illuminated,  the  Molmer  and  Sorensen 
scheme6,14  requires  that  they  all  participate  equally  in  the  inter¬ 
mediate  motional  excitation,  which  implies  that  the  only  suitable 
mode  for  arbitrary  N  is  the  centre-of-mass  mode.  However,  this 
mode  has  a  practical  disadvantage  because  in  our  experiments 
fluctuating  ambient  electric  fields  cause  it  to  heat  at  a  significant 
rate.  Although  for  large  8  the  entanglement  operation  is  largely 
independent  of  the  motion,  so  that  heating  is  unimportant,  in  the 
small-5  case  it  is  necessary  that  motional  decoherence  be  avoided. 
Modes  involving  only  relative  ion  motion  couple  to  higher 
moments  of  the  field,  so  heating  of  these  modes  is  negligible15. 
For  N  =  2  and  N  =  4,  such  modes  do  exist  in  which  each  particle 
participates  with  equal  amplitude17.  In  both  cases,  they  are  sym¬ 
metric  ‘stretch’  modes,  in  which  alternating  ions  oscillate  out  of 


Figure  1  Entanglement  scheme  for  two  particles.  Each  ion  is  initially  prepared  in  the  lj) 
internal  state,  and  the  collective  motion  of  the  pair  is  cooled  to  its  ground  state  10).  Laser 
fields  oscillating  near  o>0  +  v  and  co0  -  v  couple  the  l||)  and  Iff)  states  as  shown.  By 
detuning  the  single  transition  frequencies  by  a  small  amount  8,  the  populations  of  the  Iff  1) 
and  Ifjl)  states  are  kept  small.  Then,  by  driving  the  double  transition  for  the  appropriate 
time,  the  entangled  state  (|  ff )  -  i\  jj  ))/f2  is  created.  For  four  ions,  the  same  procedure 
generates  the  state  (|  ffff )  +  /|  ||||  We  note  that  in  the  actual  experiment,  each 
of  the  single  transitions  shown  is  itself  a  two-photon  Raman  transition,  driven  by  a  pair  of 
laser  beams;  the  entire  process  therefore  consists  of  a  four-photon  transition. 
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phase.  We  use  these  modes  here.  Excitation  of  the  centre-of-mass 
mode  still  affects  the  experiment,  as  the  motion  in  spectator  modes 
modifies  the  coupling  strength  to  the  mode  of  interest14,15.  For  this 
reason,  we  initially  cool  both  the  centre-of-mass  and  stretch  modes 
to  near  their  ground  state. 

The  experiment  was  performed  using  9Be+  ions  confined  in  a 
miniature  linear  radio-frequency  trap18,  with  the  N  ions  lying  in  a 
line  along  the  trap’s  weak  axis.  Two  spectrally  resolved  ground-state 
hyperfine  levels  comprise  the  effective  spin-half  system,  with 
IJ. )  =  \F  =  2,  mp  =  —  2)  and  |t )  =  \F  =  1,  mF  =  —  l),  where  F  is 
the  total  angular  momentum  quantum  number,  and  hmF  is  the 
projection  of  the  angular  momentum  along  the  quantization  axis 
defined  by  an  externally  applied  magnetic  field.  The  hyperfine 
splitting  coo/2tt  is  approximately  1.25  GHz.  Coherent  coupling 
between  |j)  and  ||)  is  provided  by  stimulated  Raman  transitions. 
The  two  Raman  laser  beams  have  a  wavelength  of  313  nm  with  a 
difference  frequency  near  oj0,  and  are  perpendicular,  with  their 
difference  wave-vector  lying  along  the  line  of  ions.  They  are  detuned 
~80  GHz  blue  of  the  2 Py2  excited  state,  with  intensities  giving 
■Q/27T  ~  500  kHz.  For  both  the  two-  and  four-ion  experiments,  the 
desired  stretch-mode  frequency  v/2ir  was  8.8  MHz,  giving 
»7  =  0.23 /N1/2.  The  two  driving  frequencies  required  for  the  entan¬ 
glement  operation  are  generated  by  frequency-modulating  one  of 
the  Raman  beams  using  an  electro-optic  modulator. 

After  the  entanglement  operation,  the  ions  are  probed  by  illumi¬ 
nating  them  with  a  circularly  polarized  laser  beam  tuned  to  the 
2Sin(F  =  2,  mF  =  —  2)  <-►  2Pm(F  =  3,  mF  =  —  3)  cycling  transi¬ 
tion.  Each  ion  in  |j)  fluoresces  brightly,  leading  to  the  detection  of 
~15  photons  per  ion  on  a  photomultiplier  tube  during  a  200-|jls 
detection  period.  In  contrast,  an  ion  in  ||)  remains  nearly  dark. 
Because  the  number  of  photons  detected  from  a  spin-down  ion 
fluctuates  according  to  Poisson  statistics,  in  a  single  experiment  the 
number  of  spin-down  ions  can  be  determined  with  only  a  limited 
accuracy.  For  the  data  reported,  each  experiment  was  repeated  1,000 
times  under  the  same  conditions,  and  the  resulting  photon- 
number  distribution  fitted  to  a  sum  of  poissonians  to  determine 
the  probabilities  Pj  for  j  ions  to  be  in  |j).  The  results  are  given  in 
Table  1,  and  show  that  in  both  cases,  the  probabilities  for  all  N  ions 
to  be  in  the  same  state  are  large  compared  to  the  probabilities  for  the 
other  cases.  This  is  characteristic  of  the  states  |0N),  although  the  fact 


Table  1  Characterization  of  two-ion  and  four-ion  states 


N 

Po 

Pi 

p2 

p3 

Pa 

P(tl) 

2 

0.43 

0.11 

0.46 

- 

- 

0.385 

4 

0.35 

0.10 

0.10 

0.10 

0.35 

0.215 

N  is  the  number  of  ions,  Pj  denotes  the  probability  that  j  ions  were  measured  to  be  in  |i),  and 
denotes  the  amplitude  of  the  density  matrix  element  pT Uncertainties  in  p^  and  the  N  =  2 
populations  are  ±0.01 ,  and  uncertainties  in  the  N  =  4  populations  are  ±0.02. 


that  the  middle  probabilities  are  non-zero  indicates  that  we  do  not 
generate  the  entangled  states  with  perfect  accuracy. 

In  order  to  prove  that  we  are  generating  a  reasonable  approxima¬ 
tion  to  |i //N),  it  is  necessary  to  prove  that  the  populations  of  |t-  - .  T) 
and  ||...|)  are  coherent.  In  terms  of  the  density  matrix  for  the 
system,  p,  we  must  measure  the  far  off-diagonal  element  Pf  T,i- ■  .i> 
whose  amplitude  will  be  abbreviated  pqj).  This  can  be  achieved  by 
applying  a  simple  analysis  pulse  to  the  ions  before  observing  them. 
If  the  Raman  difference  frequency  is  set  to  coo  (and  the  frequency 
modulator  turned  off),  each  ion  i  undergoes  ordinary  Rabi  oscilla¬ 
tions,  evolving  according  to  the  hamiltonian 

Hf=^(e*IT><Ul,  +  e“*li>XTI«)  (4) 

where  0  is  the  phase  of  the  difference  frequency  relative  to  that  of 
the  entanglement  pulse.  This  hamiltonian  is  applied  for  time  irl2Q 
(a  tt/2  pulse),  and  the  parity 

N 

n  =  —  iyp,  (5) 

j=0 

is  observed  while  0  is  varied.  As  seen  in  Fig.  2,  for  N ions  FI  oscillates 
as  cos  N0,  and  the  amplitude  of  this  oscillation  is  in  fact  2p(|j)  (ref. 
2).  The  resulting  values  are  given  in  Table  1.  From  the  data  shown  in 
the  table,  our  state  preparation  fidelity 

F  =  ('I'nM'Pn)  =  *  (.Pm  +  Pq)  +  Pm>  (6) 

can  be  determined,  where  Pqj  is  the  population  of  |f. .  -t)  and  Pq)  is 
the  population  of  ||...|).  For  N  =  2  we  achieve  F  =  0.83  ±  0.01, 
while  for  N  =  4,  F  =  0.57  ±  0.02. 
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Figure  2  Determination  of  pra.  a,  Interference  signal  for  two  ions;  b,  four  ions.  After  the 
entanglement  operation  of  Fig.  1 ,  an  analysis  pulse  with  relative  phase  0  is  applied  on  the 
single-ion  |1 )  <—  |  T  >  transition.  As  0  is  varied,  the  parity  of  the  N  ions  oscillates  as 


cos  (V0,  and  the  amplitude  of  the  oscillation  is  twice  the  magnitude  of  the  density-matrix 
element  Each  data  point  represents  an  average  of  1 ,000  experiments,  corresponding 
to  a  total  integration  time  of  roughly  10s  for  each  graph. 
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The  fact  that  pqj)  is  non-zero  is  still  insufficient  to  guarantee 
entanglement.  To  be  explicit,  a  system  with  density  matrix  p  exhibits 
N-particle  entanglement  only  if  no  decomposition  p  =  |  i/'jtX’Afc  I 

exists  with  all  the  {|i jj^}  factorizable  into  products  of  wavefunctions 
that  depend  on  fewer  than  N  particles.  For  example,  if  |i = 
(It ), +  11  );)/\/2  is  a  state  of  ion  i,  then  the  four-particle  state 
|i/')l  |l 11/')2|i/')3|i/')4  is  not  entangled,  but  still  has  pm)  =  1/16.  We  note 
that  these  are  the  types  of  states  studied  in  liquid-state  nuclear 
magnetic  resonance  experiments19.  Alternatively,  for  |i/>)12  = 
(lt>ilt>2  +  ll>ill>2)/V2  and  |^)34  =  (IT>3IT>4  +  U>3U>4)/V2,  the 
state  exhibits  two-particle,  but  not  four-particle  entangle¬ 

ment,  and  has  p,jj,  =  1/4. 

To  establish  that  we  are  actually  observing  N-particle  entangle¬ 
ment,  consider  an  arbitrary  factorizable  wavefunction 

h/'r)  =  [flit  •••  T )x  +  I’ll  •••  l)x  +  ••■IMt  •••  T )y  +  d\l  •••  l)y  +  ■•■] 

(7) 

where  X  and  Y  refer  to  two  distinct  subsets  of  the  N  particles,  with 
|t- . .  T)x  indicating  the  state  with  all  particles  in  subset  X  spin-up,  and 
similarly  for  the  other  terms.  Normalization  of  the  factor  wavefunc¬ 
tions  requires  |a|2  +  |fo|2  =£  1  and  |c|2  +  |d|2  1,  which  can  be 

combined  and  rewritten  as 

(Ml  -  |c|)2  +  2|«c|  +  (\b\  -  Ml)2  +  2\bd\  «  2  (8) 

Since  the  squared  terms  on  the  left  are  positive,  equation  (8)  implies 
that  |ac|  +  \bd\  1,  and  in  turn  that  (|ac|  +  \bd\)2  ^  1.  Expanding 
the  square  yields  the  desired  relation20: 

^(t)  +  ^(i)  +  2P(tu  =  2F  =£  1  (9) 

where  Pm  =  \ac\2,  Pa)  =  \bd\ 2,  and  p(tl)  =  \abcd\  are  the  previously 
defined  quantities.  Since  equation  (9)  holds  for  any  separable 
wavefunction,  it  must  also  hold  for  any  separable  density  matrix. 
Both  our  N  =  2  and  N  =  4  experiments  give  F  >  1/2,  so  the  states 
they  produce  exhibit  N-particle  entanglement. 

Quantifying  the  amount  of  entanglement  present  is  a  more 
difficult  question.  A  variety  of  measures  of  entanglement  have 
been  proposed,  but  most  are  difficult  to  calculate  even 
numerically21,22.  For  N  =  2,  there  is  an  explicit  formula  for  the 
“entanglement  of  formation”  E  as  a  function  of  p  (ref.  23).  Although 
we  have  not  reconstructed  the  entire  two-particle  density  matrix, 
the  populations  measured  place  sufficient  bounds  on  the  unmea¬ 
sured  elements  to  determine  that  E  ~  0.5.  This  indicates  that 
roughly  two  pairs  of  our  ions  would  be  required  to  carry  the 
same  quantum  information  as  a  single  perfectly  entangled  pair. 

In  the  four-ion  case,  no  explicit  formula  for  entanglement  is 
known.  The  data  do  indicate  that  our  density  matrix  can  be 
expressed  approximately  as 

p  =  0.431^X^4  1  +  0.57Pincoh  (10) 

where  |^4)  is  the  desired  state  and  p^oh  is  completely  incoherent 
(that  is,  diagonal).  These  coefficients  are  determined  directly  from 
the  value  of  pqj)  in  Table  1,  together  with  the  fact  that  no  evidence 
for  other  off-diagonal  matrix  elements  was  observed.  To  determine  a 
measure  of  entanglement,  however,  it  is  necessary  to  decompose  p  as 
a  sum  of  and  a  ‘worst-case’  factorizable  matrix  pP,  which  can  be 
accomplished  as 

p  =  0.13|ikXlM  +0.87pF.  (11) 

Note  that  equations  (10)  and  (11)  both  describe  the  same  physical 
state,  but  that  in  equation  (11),  pF  consists  of  a  specific  mixture  of 
two-  and  three-particle  entangled  states  that  is  highly  unlikely  to 
occur  in  our  experiments.  In  either  description,  it  is  clear  that  our 
state-preparation  accuracy  is  limited. 

The  source  of  decoherence  in  our  experiments  is  not  entirely 
clear,  but  evidence  suggests  that  it  is  related  to  intensity  fluctuations 
in  the  Raman  laser  beams,  a  problem  we  are  working  to  understand 


and  correct24.  The  presence  of  decoherence,  and  the  fact  that  it 
affects  the  four-ion  experiment  more  strongly  than  the  two-ion  one, 
illustrates  the  need  to  carefully  define  the  sense  in  which  our 
entanglement  operation  is  ‘scalable’.  Any  entanglement  experiment 
is  more  sensitive  to  decoherence  as  the  number  of  particles  involved 
is  increased,  unless  sufficient  accuracy  can  be  achieved  for  error- 
correction  schemes  to  be  usefully  applied.  Such  schemes  are  thought 
to  require  an  error  rate  of  the  order  of  10-4  per  operation25,  and  we 
are  certainly  far  from  this  regime.  However,  even  if  such  a  level  of 
fidelity  were  to  be  achieved,  applications  such  as  quantum  comput¬ 
ing  still  require  very  large  entangled  states  to  be  generated  in  a 
reasonable  amount  of  time  and  using  a  reasonable  amount  of 
resources.  The  method  demonstrated  here  is  important  in  this 
regard,  since  it  uses  only  a  single  operation  and  requires  a  time 
that  scales  roughly  as  N  . 

In  the  language  of  quantum  information  science,  we  have  realized 
a  four-quantum-bit  logic  gate.  This  system  is  relevant  for  the  future 
development  of  quantum  information  technology,  as  such  states 
may  be  used  to  implement  quantum  error-detection  schemes26  or  to 
make  rudimentary  demonstrations  of  quantum  algorithms27,28. 
Entanglement  of  four  particles  is  also  interesting  in  its  own  right, 
as  such  states  can  show  strong  violations  of  local  realism5.  Even  the 
two-particle  Bell’s  inequality  measurement  would  be  interesting  to 
implement,  as  the  near-perfect  detection  efficiency  for  ions  would 
eliminate  the  “fair  sampling”  hypothesis  which  has  been  required  in 
other  experiments29.  In  addition  to  improved  fidelity,  applications 
such  as  these  do  require  the  ability  to  perform  individual  manipu¬ 
lation  and  detection  of  each  ion,  but  this  is  not  expected  to  be  a 
severe  experimental  challenge30.  □ 
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Solid  hydrogen,  a  simple  system  consisting  only  of  protons  and 
electrons,  exhibits  a  variety  of  structural  phase  transitions  at  high 
pressures.  Experimental  studies1  based  on  static  compression  up 
to  about  230  GPa  revealed  three  relevant  phases  of  solid  molecular 
hydrogen:  phase  I  (high-temperature,  low-pressure  phase),  phase 
II  (low-temperature  phase)  and  phase  III  (high-pressure  phase). 
Spectroscopic  data  suggest  that  symmetry  breaking,  possibly 
related  to  orientational  ordering1’2,  accompanies  the  transition 
into  phases  II  and  III.  The  boundaries  dividing  the  three  phases 
exhibit  a  strong  isotope  effect3,  indicating  that  the  quantum- 
mechanical  properties  of  hydrogen  nuclei  are  important.  Here  we 
report  the  quantum  distributions  of  protons  in  the  three  phases  of 
solid  hydrogen,  obtained  by  a  first-principles  path-integral  mol¬ 
ecular  dynamics  method.  We  show  that  quantum  fluctuations  of 
protons  effectively  hinder  molecular  rotation — that  is,  a  quantum 
localization  occurs.  The  obtained  crystal  structures  have  entirely 
different  symmetries  from  those  predicted  by  the  conventional 
simulations  which  treat  protons  classically. 

The  structures  of  these  broken-symmetry  phases  have  been 
extensively  investigated  both  by  experimental1-3  and  by  theoreti¬ 
cal4-9  studies,  though  the  results  are  still  controversial.  In  most 
theoretical  studies,  even  when  they  precisely  compare  total  energies 
of  model  structures  on  the  basis  of  electronic-structure  calculations, 
quantum  fluctuations  of  protons  are  usually  neglected,  or  roughly 
discussed,  perhaps  due  to  some  technical  reason.  But  as  shown  by 
our  previous  study  on  impurity  muonium  (a  bound  state  comprising 


a  positive  muon  and  an  electron)  and  hydrogen  in  crystalline 
silicon10,  quantum  states  of  light  particles  may  exhibit  distributions 
significantly  different  from  those  expected  from  stability  analyses 
based  on  the  classical  potential-energy  surface.  It  is  therefore 
desirable  to  examine  quantitatively  the  quantum-mechanical 
properties  of  protons  in  solid  hydrogen,  especially  when  we  are 
concerned  with  its  symmetry  breaking  under  high  pressure.  Natoli 
et  al.u  treated  both  electrons  and  protons  with  quantum  Monte 
Carlo  methods;  their  simulations  were,  however,  performed  with 
trial  wavefunctions  based  on  the  fixed  static  configuration  of  the 
protons,  and  were  restricted  to  the  molecules  centred  on  the 
hexagonal  close  packed  (h.c.p.)  lattice  sites. 

The  first-principles  path-integral  molecular  dynamics  (FP- 
PIMD)  method  enables  us  to  incorporate  quantum-mechanical 
properties  of  protons  in  conventional  Car-Parrinello-type  first- 
principles  simulations.  The  basic  formalism  was  first  presented  by 
Marx  and  Parrinello12.  Recently  this  scheme  was  further  developed 
by  the  present  authors10,  and  we  use  this  in  the  present  study.  We 
consider  a  supercell  containing  N  (=  64)  atoms,  which  is  subject  to 
the  periodic  boundary  conditions.  To  represent  the  quantum 
properties  of  protons  in  the  path-integral  formalism,  imaginary 
time  I3h  (where  (3  is  the  inverse  temperature  in  units  of  Boltzmann 
constant)  is  divided  into  P  finite  time  slices;  each  proton  is  thus 
represented  by  a  polymer  consisting  of  P  ‘beads’  interacting  via 
intrapolymeric  harmonic  forces13.  Exchange  effects  between  pro¬ 
tons  are  neglected  in  this  study.  Assuming  the  Born-Oppenheimer 
approximation,  interatomic  forces  between  protons  at  each  atomic 
configuration  and  at  each  time  slice  are  determined  by  electronic- 
structure  calculation  based  on  the  density- functional  theory  (DFT). 
The  resultant  classical  system  of  NP  particles  is  simulated  by 
molecular  dynamics  (MD)  at  constant  volume  and  temperature 
with  the  aid  of  the  Nose-Hoover  chain  thermostat14.  In  the 
evaluation  of  the  exchange-correlation  potential  in  the  DFT  calcu¬ 
lation,  we  adopt  the  generalized  gradient  approximation  (GGA)15. 
The  interaction  between  a  proton  and  an  electron  has  been 
described  by  the  norm-conserving  non-local  pseudopotential  of 
the  Troullier-Martins  type16;  the  electronic  wavefunctions  have 
been  expanded  in  plane  waves  with  a  cut-off  energy  of  50  Ry. 
Integration  over  the  first  Brillouin  zone  has  been  achieved  by 
sampling  8  uniform  k-points.  Energy  eigenvalues  of  the  lowest  32 
occupied  bands  have  been  explicitly  calculated,  since  the  system 
should  remain  an  insulator  in  the  relevant  pressure  range1.  The 
ground-state  electron  density  has  been  obtained  by  minimizing  the 
total  energy  functional  through  the  conjugated-gradient  method17. 
We  have  verified,  through  static  DFT  calculations  with  these 
parameter  conditions,  that  the  potential  energy  surface  associated 
with  molecular  rotation  in  the  Pca21  structure  obtained  earlier  by 


Figure  1  Illustration  of  the  Pc£^  structure  projected  onto  the  x-y  plane.  Thick  arrows 
represent  molecules  on  the  a-plane  and  thin  arrows  depict  those  on  the  (3-plane,  both 
pointing  towards  the  positive-z  hemisphere.  These  two  planes  are  apart  from  each  other 
by  cl 2  in  the  z-direction.  The  quantity  6  measures  the  angle  between  a  molecular  axis  and 
the  c(z)-axis  of  the  crystal:  <p  is  the  corresponding  azimuthal  angle.  The  Cmc2i  structure 
is  obtained  by  setting  <p  =  90°. 
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